Graph Network Feature Space Fusion for Predicting Irregularly Sampled Medical Time-Series Data: Deep Learning Model Development and Validation Study

Introduction

Irregularly sampled time-series data is a common data pattern in fields such as medicine, meteorology, and transportation. Due to issues such as unequal sampling time intervals and different sampling frequencies between different features, models often use a method of sampling data at certain time intervals to align the overall data queue and fix the time interval of the sequence during the analysis and prediction of irregularly sampled time-series data [1]. For missing values in data series, existing models often use adjacent values or mean values to fill in the gaps [2]. This will cause loss of fine-grained information in time-series data and also discard the time-scale information in time-series data, which will reduce the prediction accuracy of the model when predicting and analyzing complex multiobjective downstream tasks and make it difficult to handle different types of prediction tasks simultaneously. The model constructed in this way is often difficult to solve complex and changing prediction tasks faced in the real world, especially for time-series data analysis in the medical field, which requires models with high prediction accuracy. In the medical field, it is often necessary to perform complex multitask analysis on a large amount of feature data generated in a short period of time for emergency patients, providing information for clinicians to construct treatment plans from multiple dimensions [3]. By constructing targeted treatment plans based on multiple prediction results, the mortality rate of patients can be effectively reduced.

Retrospective studies on patients’ time-series data in medicine have always been an important means to promote the development of basic medicine and explore changes in disease courses. Medical researchers can construct targeted treatment plans by discovering the changing patterns of characteristic data. The analysis methods for time-series data initially used statistical methods to analyze factor changes. With the proposal of machine learning and deep learning algorithms, researchers have built various artificial intelligence models for the classification and regression prediction of medical data [4-6]. However, due to the subjectivity of patients’ willingness to go to the hospital in reality, the sampling time interval of medical time-series data is often uneven, and there are missing values at some time points. At the same time, such data characteristics will also make the model face the problem of uneven input data. In machine learning models, researchers use the method of temporal sampling to align the data. Missing values are processed by zero-value filling or adjacent-value filling. For example, Yoon et al [7] constructed a machine learning model for time-series data of hypotension in ‌intensive care unit patients to explore the relationship between hypotension data and patient mortality. To address the issue of unequal sampling intervals, they aligned the overall data by performing mean sampling every 5 minutes. Delahanty et al [8] predicted the Sequential Organ Failure Assessment scores of patients to assess the onset time of sepsis. For missing data at certain time points, the authors filled in missing values using neighboring values. Wu et al [9] developed a machine learning model to classify the severity levels of chest pain in emergency patients. They used zero padding to handle missing values in the data. While this method introduces a significant amount of noise and homogenizes the original patient features, it allows the model to achieve reasonable accuracy for single-task applications. However, it struggles to maintain high predictive accuracy in multitask scenarios.

With the advancement of model construction methods and the introduction of deep learning algorithm frameworks, the increased flexibility in acceptable input structures has led to the development of various methods for handling missing values in medical time-series data. Lee et al [10] input each patient’s time-series data through three channels: measurement type, measurement value, and measurement behavior. They used an attention mechanism network to fuse and extract overall features, which was applied to predict patient mortality. Bardak and Tan [11] used medical order text data generated from patient clinical treatments to extract features, which were then used to supplement missing values in patients to meet the data requirements for subsequent model training. Xu et al [12] constructed a graph network data structure by applying regular sampling to patient clinical time-series data to address the impact of missing values on model training. They ultimately extracted features and fitted downstream tasks using a gated mechanism and a graph attention network. A comparative analysis of the processing methods of various deep learning models shows different strategies for handling missing values. Encoding the data into a graph network structure improves interpretability and meets the needs of medical research while minimizing the impact of missing values on the model to the greatest extent.

Aiming at the above problems, this paper proposes a data construction method based on a graph network. The feature channels of clinical time-series data of each patient and the recording time points constitute nodes in the graph network. The feature vector of the patient’s feature channel node is composed of the Kullback-Leibler (KL) divergence formed by data distribution in the feature channel, and the generated features are embedded into the feature correlation space for distribution. The patient’s clinical time-series data is input into the model for training after passing the above rule-based diagram network transformation. In the model, the point weight and edge features are extracted from the corresponding space through the model respectively, and fused under the action of a learnable head. The fused features are matched with multiobjective tasks required in clinical practice through the decoder, and the loss function is constructed according to the classification task and the regression task, respectively.

Graph network data structure is a type of data structure that can simultaneously encode spatial information and temporal change information in time-series data, which is widely used in the field of traffic and weather prediction analysis. In the early stage of research, researchers mainly focused on node features and connection relationships for classification and prediction tasks. Lu et al [13] constructed a dual-attention mechanism based on node information in traffic flow data to adaptively predict situations where nodes in the graph network change with time and space simultaneously. Jin et al [14] constructed a set of spatio-temporal prediction model of a graph network for the change of particulate matter 2.5 and predicted the change of node characteristics with time through node clustering in graph network. Sun et al [15] constructed the temperature of the ocean surface as the node feature in the graph network and constructed the transportation direction of ocean currents as the connection edge between nodes and used long short-term memory model to predict the changes of the sea surface temperature.

With the development of more model structures and computer computing power, researchers have added more information such as edge weights, edge connection directions, and point features to the graph network, which improves the final prediction accuracy of the model and can be adapted to more target downstream tasks. For example, Yin et al [16] constructed a dynamic graph network structure for traffic flow data and integrated the edge features and point features in the data into the model at the same time, so that the model could make correct predictions when the network structure changed and better fit the real changes. Feng et al [17] established the graph network of the data at each time point and also established the relationship graph network between different time points, which was used to use the delay data to guide the original data for prediction and analysis. He et al [18] proposed a multilevel fusion spatio-temporal analysis network, which fuses feature extraction from different dimensions in graph networks to adapt to downstream multiobjective tasks. Therefore, this paper uses the data structure of graph network to encode the medical clinical time-series data and makes full use of the point feature space and edge connection weights in the graph network structure to construct the features contained in the clinical time-series data. Based on the data structure of graph network, a time-series analysis and prediction model of multiobjective tasks is constructed to meet the clinical needs.

The data structure of a graph network can retain the fine-grained information and short-term feature change-in information in the patient time-series data to the greatest extent and will not bring additional noise information to the model. The construction of multiobjective tasks can enable the model to meet the needs of clinicians in the process of clinical practice, help doctors predict multidimensional changes in patients’ conditions, and formulate treatment plans for patients. The innovation points proposed in this paper can be summarized as follows:

1. This paper proposes a method to transform medical clinical time-series data into a graph network structure, which preserves the fine-grained information in the data to the greatest extent while reducing the introduction of noise information. The model avoids the problem of missing values and unequal sampling intervals in medical time-series data.
2. This paper proposes a training model for a graph network composed of time-series data, which extracts feature vectors from point feature space and edge weight space, respectively, and fuses them under the action of a learnable head to fit downstream target tasks. Both the classification task and the regression task achieve high prediction accuracy.
3. By extracting feature weights at different time points, this paper provides an interpretable analysis of critical features influencing patient survival over time, thereby assisting clinicians in making targeted and informed decisions at each stage.

Methods

Overview

In this section, we will introduce the method of constructing a graph network data structure for time-series data proposed in this paper. In the constructed graph network data, the point features and edge weights are respectively extracted and fused to construct the overall features of the patient time-series graph network. Finally, we will introduce the classifier and the overall model training framework built for the downstream multiobjective task. The overall symbols are summarized in Table 1.

Graph Network Construction

To retain the fine-grained information in the patient time-series data to the greatest extent and reduce the introduction of noise information, this paper constructs the clinical time-series data as a graph network structure to input into the model calculation. In the construction of the graph network data structure, the feature space of nodes is first constructed. By disassembling the clinical time-series data, it is found that each data measurement behavior can be divided into three dimensions: measurement time (t), measurement feature (f), and measurement value (v). Based on the calculated time range (T_m), the calculated number of feature channels (F_m), the total number of measured behaviors (V_m), and the overall time-series data of the patient can be expressed as $\left\{\right(t_1,f_1,v_1),(t_1,f_2,v_2),....(t_a,f_b,v_c\left)\right\}a\le T_m,b\le F_m,c\le V_m$. Among them, the measurement time and the measurement feature constitute the time node (P_t) and the feature node (P_f) in the graph network data, respectively, and the measurement value is constructed as the connection edge weight (E_w) between the time node and the feature node. The feature space of the feature node (P_f) is mapped in the feature correlation space and composed of the KL divergence values calculated between the original data within the computational time range across the F_m feature channels of each patient. This study references multiple studies on feature space construction based on KL divergence [19,20]. The use of KL divergence to measure distributional differences between distinct physiological variables (eg, heart rate and temperature) may indeed appear counterintuitive from a clinical perspective. However, it is important to clarify that KL divergence is not used here to infer causal or biological relationships between different features. Instead, it serves as a similarity metric to construct the feature node embeddings in the graph network. Specifically, it quantifies the distributional alignment between different physiological time series within the same patient, which may reflect the degree of systemic physiological integration or dysregulation—a key characteristic of critical illness progression. This approach is supported by prior work in multivariate time-series analysis, where distributional similarity is used to capture latent coordination among physiological systems. The specific construction method is as follows: first, select feature f; second, compute the KL divergence values between feature (f) and all other features to form the point feature vector for this feature; repeat this process until the point feature vectors for all features are constructed. The feature distributions used for a specific prediction time point are calculated strictly using data recorded strictly before that (T_m) time. For example, when constructing the graph network for prediction at time (T_m), only measurements with time stamps < T_m are used to compute the KL divergence and edge weights. No future information is incorporated, thereby preventing temporal data leakage. The feature space integrates the information correlation and distribution correlation between feature channels and establishes the relationship between the edge weight value and the point feature space. At the same time, the KL divergence values in the overall feature space are normalized to ensure that the feature distribution is in the feature correlation space and reduce the impact on model training. The calculation method of the feature space (P_f) is as follows:

$S_{pf}=\{F_1^{pf},F_2^{pf}......F_b^{pf}\}b\le F_m$(1)

${\displaystyle {\displaystyle F_b^{pf}=\left[KL\left(f_b,f_1\right),KL\left(f_b,f_2\right).......KL\left(f_b,f_{F_m}\right)\right]b\le F_m}}$(2)

Probability distribution construction for KL divergence under sparse data: for each feature channel, let its available measurements before prediction time be $\{v_1,v_2,...v_n\}$. This paper constructs a discrete probability distribution over a fixed set of M bins spanning the global range of all features. When n=2, we use kernel density estimation with a Gaussian kernel; when n=1, the distribution is a Gaussian centered at that single value; and when n=0, we assign a uniform distribution. The KL divergence between 2 features is then computed on these discrete distributions. Where, S_pf represents the overall spatial feature set, and F^pf_b represents the feature vector of a feature node. Meanwhile, the feature space mapping of the time mode P_t is in the temporal difference embedding space, built by the difference between the time node and the measured time of the remaining nodes. The feature space is used to reflect the correlation of the time values in the measured time series. The calculation time point is strictly limited within the calculation time range. The information provided by the time node can assist the model in mining the time scale between different measurements and provide the change information of the time scale for the model training. The calculation method of the feature space P_t is as follows:

${\displaystyle {\displaystyle {\mathrm{S}}_{pt}=\left\{{\mathrm{F}}_1^{pt},{\mathrm{F}}_2^{pt},\dots ,{\mathrm{F}}_a^{pt}\right\},a\le {\mathrm{T}}_m}}$ (3)

${\displaystyle {\displaystyle \begin{array}{r}{\mathrm{F}}_a^{pt}=[-\Vert ({\mathrm{t}}_a,{\mathrm{t}}_1)\Vert ,-\Vert ({\mathrm{t}}_a,{\mathrm{t}}_2)\Vert ,\dots ,\\ -\Vert ({\mathrm{t}}_a,{\mathrm{t}}_m)\Vert ],a\le {\mathrm{T}}_m\end{array}}}$ (4)

Where, S_pt represents the overall spatial feature set, and F^pt_a represents the feature vector of a time feature node. The connection edge E_w of the graph network is composed of v_c in the measurement behavior, and the existence of the connection edge is directly related to the existence of the measurement behavior. In this way, the fine-grained information in the data is retained to the greatest extent, and the introduction of noise information is eliminated. The construction method of the graph network is shown in Figure 1.

‎

Model Construction

After constructing a graph network with time-series data, this paper proposes the construction of a graph convolutional network (GCN) by using the point feature space and the edge weight connection relationship to extract the overall features in the point feature space. At the same time, a multihead attention mechanism network with a gated mechanism connecting multiple layers is proposed for feature extraction in the edge weight feature space. Finally, a learnable weight head is added to the features of the 2 dimensions for feature fusion to meet the downstream multiobjective tasks.

First, in the GCN, the feature spaces S_pf and S_pt are jointly constructed as the feature matrix X_f of the nodes. The collected values of each feature for an individual patient at different time points are normalized, and the connection values corresponding to time nodes and feature nodes form the edge weight matrix A_w, and the learnable parameter Wis added. The GCN [21] network with L) layer input edge weights is used to construct and calculate the overall feature F_p of the node feature space. The calculation is as follows:

$H\left(1\right)=\sigma \left(A_w{X}_f{W}_0\right)$(5)

${\displaystyle {\displaystyle F_p=H\left(L\right)=\sigma \left(A_{w}H\left(L-1\right)W_{L-1}\right)}}$(6)

Second, during the computation process of the GCN network, the model primarily focuses on the intrapatient variations, and the edge features tend to gradually weaken as the training depth increases. To fully express the features in the edge weight space, this paper proposes a targeted approach. First, each feature of the edge weight (E_w) is independently normalized across the entire set of patient samples. This normalized edge weight data is then used to construct the feature value matrix (W_E). Finally, (W_E) is input into the convolutional gated attention mechanism network.

During the computation process of the convolutional gated attention mechanism network, first, W_E is fed into the forgetting gate to construct the partially forgotten feature matrix (W_EF). Then, W_EF and W_E are input into the 2D convolution layer at the same time, and the multidimensional data (C_EF) and (C_E) are generated after convolution and input into the multihead attention mechanism layer [22] to obtain A_E(1). It is constructed as shown in the following formula:

${\displaystyle {\displaystyle W_{EF}=Sigmoid\left(h_1{W}_E+b_1\right)}}$(7)

${\displaystyle {\displaystyle C_{EF}=Conv2d\left(W_{EF}\right)}}$(8)

$C_E=Conv2d\left(W_E\right)$(9)

$A_E\left(1\right)=Attention(C_{EF},C_E,C_E)$(10)

Here, (h₁) represents the weight matrix of the first layer and (b₁) represents the bias term of the first layer. The feature A_E(1) obtained from the first layer is reinput into the convolutional gated network to generate the C_AE after passing through the average pooling layer. C_AE, A_E(1), and C_E are again input into the multihead attention mechanism layer. Repeat the operation of feature extraction through (L) layers, and finally obtain the extracted feature weight (F_e) of the edge feature space. The overall calculation formula is as follows:

${\displaystyle {\displaystyle \begin{array}{r}{\mathrm{C}}_{\mathrm{A}\mathrm{E}}(\mathrm{L})=\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{v}2\mathrm{d}(\mathrm{S}\mathrm{i}\mathrm{g}\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{d}({\mathrm{h}}_L\mathrm{A}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}({\mathrm{A}}_E(\mathrm{L}-1))\\ +{\mathrm{b}}_L))\end{array}}}$(11)

${\displaystyle {\displaystyle A_E\left(L\right)=Attention\left(C_{AE}\left(L\right),A_E\left(L-1\right),C_E\right)}}$(12)

${\displaystyle {\displaystyle F_e=AveragePooling\left(A_E\left(L\right)\right)}}$(13)

Finally, the feature (F_p) obtained from the point feature space and the feature (F_e) obtained from the edge weight space were added to the learnable weight head (H_fp) and (H_fe), respectively, to form the overall feature extracted from the time-series data graph network space. The features are computed as follows:

$F_{total}=H_{fp}{F}_p+H_{fe}{F}_e$(14)

Model Training

After obtaining the overall features extracted from the temporal graph network, we build different decoders to satisfy the downstream multiobjective task. For the classification prediction task, it is necessary to trace the overall feature changes of patients and conduct a comprehensive analysis. For the 2 different classification tasks of patient mortality and patient clinical representation, the overall feature matrix (F_total) is decomposed and arranged by time and mapped to the corresponding dimension through the feedforward neural network for training. The binary cross-entropy loss function is used to calculate the task loss value, and the calculation formula is as follows:

${\displaystyle {\displaystyle {\displaystyle \begin{array}{rl}\mathrm{B}\mathrm{c}\mathrm{e}\mathrm{L}\mathrm{o}\mathrm{s}\mathrm{s}& =-\frac{1}{N}\sum _{i=0}^N{\mathrm{y}}_i\log (\mathrm{p}({\mathrm{y}}_i))\\ & +(1-{\mathrm{y}}_i)\log (1-{\mathrm{y}}_i)\end{array}}}}$(15)

For the regression task, the key consideration is the feature vector of patients at the final time point. For regression prediction tasks such as predicting the future length of hospital stay or the future change of different physical sign data, this paper uses the features of the last dimension of the point feature space and edge weight space, which are weighted by H_fp and H_fe and then mapped to the corresponding dimension by the feedforward neural network for training. The calculation formula is as follows:

$L1loss=\frac{1}{N}{\sum }_{i=0}^N\left|y_i-f\left(y_i\right)\right|$(16)

The length of hospital stay was normalized before model training using minimum-maximum scaling. Specifically, let (y_st) denote the original stay time in hours. The normalized value (y_norm) is computed as:

$y_{norm}=\frac{y_{st}-y_{min}}{y_{max}-y_{min}}$

The overall model construction and training are shown in Figures 2A and 2B.

‎

Ethical Considerations

This research was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of the First Affiliated Hospital of Wenzhou Medical University (approval KY2023-R061). The requirement for written informed consent was waived by the Institutional Review Board due to the retrospective nature of the study. The study data are anonymous or deidentified. Meanwhile, no identification of individual images in the paper or additional material is possible.

Results

Dataset Introduction

In the experimental data part of this paper, the time-series data of patients extracted from the public datasets Medical Information Mart for Intensive Care (MIMIC) III [23,24] and MIMIC-IV [25] are used. Under the guidance of the clinical practice experience of clinicians, this paper selects 42 important features of concern in clinical practice for model construction. Meanwhile, the patient’s survival, length of hospital stay, and clinical phenotypes are constructed as the outcome variables of the model. Since the recording time of the MIMIC-IV dataset is in the transition period of characterizing the coding rules, it is split into two parts: ICD-9 (International Classification of Diseases, 9th Revision) and ICD-10 (International Classification of Diseases, 10th Revision). The data of patients collected by the Critical Care Medicine Department of the First Affiliated Hospital of Wenzhou Medical University from 2012 to 2021 are also used to verify the effectiveness of the model. The data of 2000 patients in each dataset are selected in a 1:1 ratio of death to survival to explore the effectiveness of the proposed model in a balanced and fair way. Overall statistics for the data are provided in Table 2. Among them, the format of categorical data is expressed as “number of patients to number of nonpatients,” and the format of numerical data is expressed as “mean (SD).” To verify the stability and effectiveness of the model, the proposed model was evaluated on an expanded dataset of 4000 samples with a 1:1 ratio of death to survival, analyzing temporal feature data more than 6-hour and 12-hour time windows. The results were compared with those of state-of-the-art methods, and the detailed comparison is provided in Multimedia Appendix 1.

Classification Prediction Experiment

In this paper, the proposed model is compared with advanced time-series prediction models in two single- and multi-objective tasks: predicting the final death of a patient and predicting whether multiple phenotypes of a patient occur. The study explores the advantages and scope of application of the proposed model. The experimental equipment adopted in this paper includes an Intel(R) Xeon(R) Silver 4210R CPU running at 2.40 GHz, an NVIDIA GeForce RTX 3090 graphics card with 24 GB of video memory, and 16 GB of system memory. The proposed model in this paper is compared with state-of-the-art time-series prediction models designed for irregularly sampled data, including time-aware long short-term memory [26], gated recurrent unit with decay [27], multitime attention networks for irregularly sampled time series [28], and ContiFormer [29]. Additionally, we adopt time-series analysis methods such as fully convolutional networks [30], temporal convolutional networks [31], neural basis expansion analysis for interpretable time series forecasting [32], Crossformer [33], and DSformer [34], which are trained on regularly sampled data with missing values imputed by the mean. The final comparison results provide a more multidimensional perspective for evaluating time-series prediction performance. For the time range of data length, this paper refers to the experience in the clinical practice process to construct 5 important time period lengths for training and comparison. The hyperparameter tuning of the model refers to multiple state-of-the-art hyperparameter optimization methods [35-39], and the optimal training parameters for each model are obtained. Detailed data are provided in Multimedia Appendix 2. The accuracy and area under the curve for predicting death for each model are provided in Tables 3 and 4.

From the results, it can be seen that the proposed model is in a leading position in most tasks of the 4 datasets, and the results of the proposed method are equal to or slightly inferior to the results of the comparison methods in some time periods of some datasets. Meanwhile, the proposed model demonstrates accuracy in predicting mortality on a larger-scale dataset. The detailed comparison results are provided in Multimedia Appendix 1. From the analysis, it can be found that death is a simple binary classification task, and most of the time-series analysis methods can effectively capture the change law of time-series data to make the correct judgment. The wrong part of the judgment may be due to changes in the patient’s condition, resulting in the patient’s condition getting better or suddenly worse. Moreover, in the public datasets, the clinical time-series data of patients are collected with high frequency and high data density, and the results of various time-series prediction methods are basically similar. However, in the private dataset, due to different clinical policies, the data collection frequency of patients is relatively sparse, and there are many missing values in the normal time-series arrangement, which leads to poor and unstable prediction results. The proposed algorithm can predict stably and efficiently even when there are many vacant values.

For the multitarget classification task of predicting whether a patient‘s phenotype occurs, the average accuracy, area under the curve, and precision-recall area under the curve (PRAUC) of each target are used as evaluation metrics. PRAUC is particularly suitable for imbalanced classification tasks, as it focuses on the model’s performance on the minority class—a common scenario in clinical phenotype prediction where positive cases are often rare. The prediction performance of each model is summarized in Tables 5-7.

From the Tables 5 and 6, it can be seen that the proposed model has achieved superior results in different datasets and different time periods, which fully reflects the advantages of the proposed algorithm in the face of complex downstream tasks. As shown in Table 7, the proposed method consistently achieves the highest PRAUC across most datasets and time windows, indicating its robustness in handling severe class imbalance. Meanwhile, the proposed model demonstrates accuracy in predicting multiple clinical phenotypes on a larger-scale dataset. The detailed comparison results are provided in Multimedia Appendix 1. Comparing different time-series analysis algorithms, the sampling-and-filling method does cause the loss of fine-grained information. However, the data input structure of the graph network proposed in this paper can meet the required input of the model while retaining the fine-grained information in the time-series data to the greatest extent. At the same time, the experimental results show that the simpler the model structure may produce better results when computing time-series data with sampling and filling.

Regression Prediction Experiment

In the experimental architecture of this paper, predicting the length of hospital stay of a patient is set as a regression prediction analysis task to explore the prediction accuracy of the model in the face of different types of downstream tasks. At the same time, the model can provide multidimensional information required by clinical practice. In this paper, the normalized length of hospital stay of patients is set as the prediction target, and the mean absolute error (MAE) value between the true value and the predicted value is used as the standard to evaluate the quality of the model. The comparison of regression predictions is shown in Table 8. For each dataset, the corresponding minimum and maximum values were used for scaling. Detailed dataset-specific normalization parameters are provided in Multimedia Appendix 3.

From the analysis of the table, it can be found that the prediction accuracy of the model proposed in this paper is in the leading position in most datasets and time periods and is slightly inferior to the best prediction results in a small number of time periods. This finding fully reflects the stability of the model proposed in this paper in the face of different downstream tasks. Meanwhile, the proposed model demonstrates accuracy in predicting the length of hospital stay on a larger-scale dataset. The detailed comparison results are provided in Multimedia Appendix 1. Compared with other models, it can be found that some time-series analysis models have good accuracy in classification experiments, but have large deviation values in regression prediction tasks. In this paper, the prediction indicators of different time periods of each dataset are averaged, and the regression prediction and classification prediction are projected into the same set of score space to verify the stability of the model in the face of multitask learning.

Ablation Experiment

In the ablation experiment of the model, the point feature space and edge weight space proposed in this paper are experimented on using the above 3 regression and prediction analysis tasks (death, phenotype, and stay time) on 4 datasets with different time lengths, and the contribution of feature vectors extracted from different feature spaces to the overall model is explored. Meanwhile, this paper also analyzed the connection results of 2 feature spaces under the learnable head and the direct connection results of 2 feature spaces. Results from top to bottom, the accuracy for predicting death, the average accuracy for predicting the occurrence of different representations of patients, and the MAE for predicting the length of hospital stay were arranged. The values are given in Table 9.

From the analysis of the table, it can be found that both the point feature space and the edge weight space can achieve certain prediction accuracy in the classification task. However, due to the characteristics of the graph network space and the fact that the GCN network is easy to fall into local graph space, the features extracted from the point feature space lead to unstable prediction results. The lack of measurement time information in the edge weight space also leads to low prediction accuracy of the extracted features. The combination of the two can effectively construct a stable and high-accuracy prediction model. In regression prediction analysis, the information provided by the point feature space is limited, and the results of regression prediction have large deviations. Under the action of the learnable head, it can effectively fuse with the feature vector extracted from the edge weight space. It improves the prediction accuracy and prediction stability of the model. Meanwhile, a single-head attention mechanism network and a single-layer GCN are supplemented to verify the effectiveness of hyperparameter configuration and the construction of the multilevel network framework. The results are shown in Table 10.

Discussion

Principal Findings

Aiming at the problem of uneven sampling time intervals and a large number of missing values in medical time-series data, this paper proposes to transform medical time-series data into a graph network structure for input. The feature vectors are extracted by constructing algorithms respectively for the point feature space and the edge weight space in the graph network data structure and are fused to form the feature vector of the graph network space under the action of the learnable head. The fused features are then decoded by the decoder to cope with different downstream tasks. In the exploration of model interpretability, this study extracts each feature at different time points from the trained model. The extracted features are fitted to different target tasks using logistic regression, and Shapley Additive Explanations analysis is used to obtain the changes in the contribution of features to the final target task across various time points. Figure 3 illustrates the variation of feature contributions to mortality prediction at different time scales based on 48-hour data from the MIMIC-IV ICD-10 dataset. Results for the entire dataset are provided in Multimedia Appendix 4.

‎

As analyzed from Figure 3, heart rate and glucose indicators are important factors affecting patient survival in the early stage of hospitalization. With the passage of time, carbohydrate metabolites and patient age gradually become the dominant factors.

Comparison to Prior Work

In this paper, the proposed model is tested in 5 different data time ranges in 4 datasets, and it is found that the proposed model has good prediction accuracy and prediction stability in classification and regression prediction tasks. Meanwhile, paired DeLong tests were performed on the mortality probability prediction results of the proposed model and comparative models in the MIMIC-IV ICD-9 and private datasets. The results showed that the P value of the test in MIMIC-IV ICD-9 was <0.001 and that in private data was <0.05. Additionally, paired U-tests were conducted on the MAE values of the regression predictions for length of hospital stay, revealing that P values in both MIMIC-IV ICD-9 and private data were <0.05. The specific value of P is provided in Multimedia Appendix 5. These findings fully demonstrate the effectiveness of the proposed model in improving prediction performance. Moreover, this paper conducts ablation experiments on the point feature space model and the edge weight space model and finds that the accuracy and stability of the classification prediction of the model come from the simultaneous contributions of the 2 spaces. The accuracy of regression prediction comes from the adjustment and fusion of the 2 spatial feature vectors by the learnable head.

In comparing our proposed model with existing methods, it is important to acknowledge a confounding factor in the baseline experimental design. Our graph-based model operates directly on raw, irregularly sampled time-series data without any imputation, whereas standard deep learning models, such as a fully convolutional network and temporal convolutional network, are trained on data that have been preprocessed with mean imputation to enforce regular sampling intervals. This design choice conflates 2 distinct sources of potential performance gain: the architectural advantage of the graph network and the representational benefit of avoiding information loss due to imputation. The primary rationale for this comparison was to reflect a realistic clinical deployment scenario, where missing values are prevalent, and imputation remains the default preprocessing step for conventional models. As such, the current benchmark captures the end-to-end practical advantage of a model that inherently accommodates irregular data without the need for imputation. However, we acknowledge that this setup does not isolate the contribution of the graph architecture alone. Despite this limitation, we believe the current comparison remains clinically meaningful, as it reflects the performance gap between models that require imputation and those that do not in real-world settings.

Limitations

The model proposed in this paper also has the following limitations: (1) since the method proposed in this study fully uses the information in the data to construct nodes and edges and uses both nodes and edges for computation simultaneously, it imposes certain requirements on the training space; and (2) as the number of computational features and the scope of computation time increase, the demand for computational space grows exponentially.

Future Directions

In the future, this study intends to develop a graph network construction paradigm with a more compact structure, which can be compatible with longer time-series data and more features, thereby enabling the tracking and monitoring of patients’ overall disease progression. While the proposed model demonstrates strong predictive performance on irregularly sampled medical time-series data, its scalability to large-scale clinical databases remains a challenge. As noted in the Limitations subheading, the current graph construction and joint computation over node and edge feature spaces impose significant memory overhead that grows exponentially with the number of features and time points. To facilitate real-world deployment in hospital informatics environments that routinely handle tens of thousands of patients, several optimization strategies can be explored. First, sparse graph representations could be adopted to retain only clinically meaningful edges, reducing the memory footprint. Second, feature or temporal down-sampling techniques may help control the node count while preserving essential information. Third, distributed computing frameworks or graph mini-batching strategies could partition the computation across multiple devices, enabling parallel processing of large patient cohorts. Finally, memory-efficient graph convolution operators, such as those based on neighborhood sampling or low-rank approximations, offer promising avenues for scaling graph networks to larger datasets. In future work, we aim to investigate these directions to enhance the scalability and clinical utility of the model without compromising predictive accuracy.

Conclusions

Based on the model’s prediction results, analyzing patients’ clinical time-series data from the early time period can reveal the disease progression for many patients. In addition to carefully examining the characteristic representations of patients identified by the model as having higher mortality risk, it is also possible to focus on predicting treatment regimens that are likely to be misapplied to a given patient. Through the similarity of the patient’s graph network to form the edge and through the hierarchical clustering method, a good treatment plan is mined, and the wrong treatment plan is excluded. Overall, this paper proposes a deep learning–based temporal data analysis model for accurately predicting the survival outcomes of critically ill patients, while also providing critical features that require close attention at various time points.